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Abstract. We present a novel method to compute the phase space distribution in the 
nonequilibrium stationary state of a wide class of mean-field systems involving rotators 
subject to quenched disordered external drive and dissipation. The method involves a 
series expansion of the stationary distribution in inverse of the damping coefficient; the 
expansion coefficients satisfy recursion relations whose solution requires computing a matrix 
where about three quarters of the elements vanish, making numerical evaluation simple and 
efficient. We illustrate our method for the paradigmatic Kuramoto model of spontaneous 
collective synchronization and for its two mode generalization, in presence of noise and 
inertia, and demonstrate an excellent agreement between simulations and theory for the 
phase space distribution. 
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1. Introduction 

Spontaneous collective synchronization in a large population of coupled oscillators of varying 
frequencies occurs in a wide variety of systems spanning length and time scales of several 
orders of magnitude. Examples are yeast cell suspensions [T], cardiac cells [2], hreflies [3], 
Josephson junctions [1], atomic recoil lasers [3], animal flocks [3], pedestrian motion on 
footbridges [7] , audience applause in concert halls [8] , and many others [9] . The paradigmatic 
minimal model to study synchrony and its emergence from asynchronous/incoherent phase 
is the celebrated Kuramoto model, involving phase-only oscillators of distributed natural 
frequencies coupled via a mean held mm- The model enjoys the status of one of the 
most studied nonlinear dynamical systems, which continues to provoke many unresolved 
issues [T2] . 

In the Kuramoto model, the oscillator phases follow a hrst-order evolution in 
time. A generalized second-order stochastic dynamics, initially studied to model better 
synchronization among hashing hrehies, accounts for the hnite moments of inertia of 
the oscillators (thus, the oscillators become instead the rotators), and for the stochastic 
huctuations of the natural frequencies in time [T3HT5]. The generalized dynamics without 
stochasticity also arises in electrical power distribution networks [MlIIT]. 

In a diherent context, one may interpret the generalized Kuramoto dynamics as that 
of a system of interacting particles driven by quenched disordered external torques, and 
evolving in presence of an external heat bath and dissipation; the Kuramoto model is 
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recovered as the overdamped, noiseless dynamics. This interpretation offers the very 
promising possibility of studying the dynamics from statistical mechanical perspectives 
by using tools of kinetic theory [181 [19]. The generalized model is a rich laboratory to 
study many-body nonequilibrium dynamics in presence of external drive and quenched 
disorder. The system relaxes at long times to a nonequilibrium stationary state (NESS) 
[mils]; the synchronized state (respectively, the incoherent state) corresponds to a spatially 
inhomogeneous (respectively, spatially homogeneous) NESS. 

Unlike in equilibrium, a NESS is characterized by a violation of detailed balance, and 
an associated phase space distribution that cannot be expressed in the Gibbs-Boltzmann 
form ~ exp{—(3H) in terms of a Hamiltonian H and an inverse temperature [3. Instead, the 
NESS distribution has to be obtained by studying on a case-by-case basis the underlying 
dynamical model. This task proves daunting especially for many-body interacting systems; 
exact results are known only for simple models, often via tour de force |2Q|, while for more 
complex models, simulations and approximation methods are invoked IZH. 

In this paper, we provide a novel method to obtain the inhomogeneous NESS single¬ 
rotator phase space distribution for the generalized Kuramoto model in the thermodynamic 
limit. The resulting distribution has a nontrivial form with respect to the Gibbs-Boltzmann 
distribution in equilibrium. Our result provides an example of a computation of a nontrivial 
probability distribution associated to a NESS in presence of quenched disorder. Significantly, 
our proposed method applies not just to the case at hand, but, in fact, to any periodic two- 
body mean-field interaction potential, thereby providing a general framework to compute 
NESS distribution in a wide class of systems. 

Our method involves a series expansion of the stationary distribution in terms of the 
inverse of the damping coefficient, as opposed to an expansion in the same parameter of 
the Kramers operator for the time evolution of the single-rotator distribution [22]. In 
our method, the expansion coefficients satisfy recursion relations whose solution requires 
evaluation of a matrix where about three quarters of the elements vanish, making the 
numerical implementation very efficient computationally and also simple compared to the 
operator-expansion method. Also, in contrast to application of the latter method to systems 
in external fields m, we develop our method for mean-field systems that require self- 
consistent evaluation of the mean fields. Note that the series expansion of the distribution 
is asymptotic in nature (just as the expansion of the Kramers operator [22]), and one can 
resort to the Borel summation method to sum the series properly [23]. In the context of the 
generalized Kuramoto dynamics, we demonstrate an excellent agreement of our theory with 
simulation results for the phase space distribution. The latter is characterized by a spatially 
non-uniform temperature profile distinctive of NESSs and absent in equilibrium. 

We remark that in a recent work [18] we studied the generalized Kuramoto dynamics, 
with a focus on obtaining the complete phase diagram of the model for a general unimodal 
frequency distribution of the oscillators. Based on simulations and some analytic bounds. 
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we demonstrated that the model shows a nonequilibrium first-order phase transition from 
a synchronized phase at low values of the relevant parameters to an incoherent phase at 
high values. However, our analytical approach allowed to obtain only the homogeneous 
NESS single-rotator phase space distribution and to study its dynamical stability, while we 
could not compute analytically the distribution of the inhomogeneous NESS, namely, the 
synchronized state. In the present work we do not address the issue of phase transitions, but 
we focus on a method that computes the inhomogeneous NESS, not just for the Kuramoto 
potential, but for any other periodic mean-field interaction potential. Thus, we also fill the 
gap left in our previous work. 

The paper is organized as follows. In the following section, we state the setting of 
the class of mean-field models that forms the object of our study, and write the dynamical 
equations in a convenient dimensionless form. In Section |3l we discuss the Kramers equation 
for the time evolution of the single-rotator phase space distribution, and in particular, 
present in detail our method to compute its stationary solution for the inhomogeneous phase. 
Some of the technical details are relegated to the Appendix. In Section 01 we illustrate our 
method by considering the representative case of the Kuramoto interaction potential, and 
demonstrate an excellent agreement between theory and simulations for several physically 
relevant observables. The paper ends with conclusions. 


2. The class of models 


We now turn to deriving our results, by hrst stating the setting of the generalized Kuramoto 
dynamics: We have N globally coupled rotators of same moment of inertia m, with the Ah 
rotator, i = 1,2,..., N, having its natural frequency Ui a quenched random variable given 
by a common distribution Q{u). The phases 6 *j’s G [0,27r] and the angular velocities Uj’s 
follow the equations [TdlfTK] 


dt 


= Vi 


( 1 ) 


N 


dvi K duiO 

m— = --fVi + ^ 2^ - 


j 


dt 


i=i 


d9i 


VDT]i{t), 


where 7 is the damping coefficient, K is the coupling constant that is scaled down by N to 
have a well-defined behavior of the associated term in the thermodynamic limit N ^ 00 , 
u{6i — 9j) is the two-body mean-held interaction potential [23], while 7 ^ is a Gaussian white 
noise: 


(hiW) = 0 , = 26ij6{t - t'). ( 2 ) 

Here, the angular brackets denote averaging with respect to noise realizations. The constant 
D in Eq. ([T]) quantihes the strength of the noise force. Noting that u{9) is periodic and even 
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( 3 ) 


in 9 [25], and taking m( 0) = 0 without loss of generality, a Fourier expansion yields 

CX) 

u{6) = — cos(s6')]. 

5=1 


The Kuramoto potential corresponds to the choice Ui = 1,Ms>i = 0. The noiseless {D = 0), 
overdamped (m/y —?■ 0) limit recovers the Kuramoto dynamics [TUlfTT] : 

N 

0^), (4) 


dOi K ^ 


i=i 


where K = iF/y. 

Common to studies of the Kuramoto model, we consider a unimodal Q{oj), i.e., one 
which is symmetric about a single maximum (same as the average (oj)), and with width 
a. The dynamics ([1]) is invariant under 9i ^ 9i + {u)t,Vi Vi + Ui + {oj), and 

the effects of a may be made explicit by replacing Ui in the second equation with era;*. We 
consider from now on the dynamics ([I]) with cu* —)■ aui, and take to have zero mean 
and unit width, without loss of generality. 

On interpreting the model ([I]) as that of interacting particles, m becomes the mass, 9i 
the angular coordinate for the motion along a unit circle, Vi the angular velocity, and yo;* 
the quenched disordered external torque [TSlIT^ . Introduced to mimic fluctuations of the 
natural frequencies [26|, the Gaussian noise can also be interpreted as fluctuations due to 
coupling to a heat bath at temperature T, and one may invoke the fluctuation-dissipation 
relation to relate the strength of the noise to the temperature T, a.s D = yfc^T [27|. In 
this case, Eq. ([1]) in the absence of the uj^s describes the dynamics of a Hamiltonian system 
in contact with a heat bath, and the stationary state is in equilibrium, with probability of 
conhgurations ~ exp{—H/T), where H is the Hamiltonian 

N 

9,), (5) 


H = V ^ + — V u{9i 
^ 2m 

i=l 


with Pi = mvi the angular momentum of the Tth particle. In presence of cuds, the dynamics 
relaxes to a NESS [T8|[T9]. 

It is convenient for further analysis to write Eq. ([1]) in a dimensionless form. Introducing 
dimensionless variables 


t = t\/ K/m, 

(6) 

Vi = Vi^/m/K, 

(7) 

l/V^ = ■y/V Km, 

(8) 

a = ya! K, 

(9) 

T = T/K, 

(10) 

Viit) = 

(11) 
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Eq. ([T]) becomes 


dt 


= Vi, 


where 


dVj 

dt 



OO 

+ sUgRs sin(s-^ — s9i) -\- aui + 

S = 1 



( 12 ) 


3 = ^ 

(13) 

and the noise satishes 


{Viit)Vjit')) = 26 ij 6 {t-t'). 

(14) 

For m = 0, using dimensionless time 


t = t{K/-f), 

(15) 

the dynamics ([T]) becomes the overdamped motion 


(^0. , — 

sUgRs sin(s'0 — s9i) + aoji -|- yTri^(t), 

L16 

(16) 


which in the case = 1, Ms>i = 0 and at T = 0 becomes the Knramoto dynamics. Associated 
with the sth Fonrier mode of the interaction potential is the magnitnde of the mean held Rs 
acting on one rotator dne to its interaction with all the others, while s'ljj is the corresponding 
phase. In particular, Ri measures complete phase coherence among all the rotators, and its 
stationary value Rf = Ri{t ^ oo) serves as the synchronization order parameter. From now 
on, we consider the dynamics flT^ that, besides N, depends on the parameters m, T, and 
a; we will drop the overbars for notational simplicity. Note that dynamics (lT2l) reduces to 
dynamics flT^ in the overdamped limit. 

In the next section, we discuss the Kramers equation for the time evolution of the single¬ 
rotator phase space distribution, and in particular, our method to compute its stationary 
solutions. 


3. The Kramers equation and its stationary solutions 

In the thermodynamic limit, the dynamics flT^ is described by the single-rotator phase space 
distribution f{6,v,uj,t), giving at time t for each u the distribution probability of rotators 
with phase 6 and angular velocity v. We have 

f {6,V,uj,t) = f {6+ 271,v,uj,t), (17) 
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and the normalization 

J d9dv f {9, v,u),t) = 1 V cj. (18) 

The time evolntion of /, obtained by trnncating to lowest order in 1/N the Bogolinbov- 
Born-Green-Kirkwood-Yvon (BBGKY) hierarchy equations, follows the Kramers equation 

ra 

df _ _ af T (fj a 

dt ^ 89 y/m dv"^ dv 
where 

= / d9dvdu g{uj)e^^^f(9,v,u,t). (20) 


V 

m 


sUgRs sin(s'0 — s9) — au \ f 


s=i 


(19) 


We are interested in the stationary state solutions of the Kramers equation, obtained by 
setting the left hand side of Eq. (IT^ to 0. As already mentioned, the stationary state is a 
NESS, unless a = 0. In the stationary state Rg and ^ are time independent; ^ can be set 
equal to 0 by redehning the origin of the 9ds, while with Rf we will denote the stationary 
value of Rg. We will also use the dehnition 

CX) 

G{9) = sUsRf sm{s9). (21) 

s=i 


The 6'-independent solution characterizing the incoherent phase, for which Rf = 0 V s and 
thus G{9) = 0, is given by [14] : 


n{9,v,u) 





( 22 ) 


The synchronized phase distribution for a = 0 is given by the Gibbs-Boltzmann measure 
~ exp[—n^/(2T) — f d6'G(0)]; for general a, we expand it as 


(23) 


where the functions bnS satisfy bn{9, cj) = bn{9+27r, cu), while $„(aa;) is the Hermite function: 

^Jax) = J ^ exp Hn{ax), (24) 

z^nlx/Ti z 


with Hn{x)^s being the n-th degree Hermite polynomial. The functions <l>„’s are orthonormal: 
J dx <hm(aa;)<h„(ax) = 6mn- Normalization of /®y°(6*, n, a;) implies that fQ"" d9 bQ(9,uj) = 1, 
while the self-consistent values of the parameters Rf are given by 


f27r 


Rf = du g(u}) / d9 bo{9,u)e 


is6 


(26) 
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Furthermore, using f da; a;<Fo(o^)®n(a3;) = l/(\^a)Sn,i, we obtain that 
[ du V{9, V, ui) = VTbi{9,u)). 


( 26 ) 


On the other hand, integrating over v the stationary state Kramers equation, we obtain that 
J dn n, cj) and, hence, bi{9,u)), does not depend on 9. 

The choice of the Hermite functions in the expansion fl2^ is motivated by the fact for 
a = 0, the distribution fsyn{9,v,u)) should have the Gibbs-Boltzmann form, fsyn{9,v,u) ~ 
exp[—(2T) — J d6'G(6')]. As will be shown later, the expansion coefficients bn for this case 
satisfy bQ{9, 0) ~ exp[— J d6'G(6')], bn{9, 0) = 0 for n > 0, so that only the n = 0 term in the 
expansion (12^ needs to be taken into acconnt; then, with <Fo(a;) ~ exp(—x^/2), the prodnct 
( vir) (v^) expansion correctly reprodnces the velocity-part of the 

distribntion ~ exp[—n^/(2T)]. 

Plugging the expansion fl2^ into the stationary state Kramers eqnation, using the known 
recnrsion relations for the Hermite polynomials, and eqnating to zero the coefficient of each 
we get 

/-Tf;dbn-l{9, u) r ——dbn^\{9,Uj) 

^ BB + ae + 


77 / 77 

bn[9,bj) W —bn-\i9,bS)\G[9) - crw] = 0 
m VI 


(27) 


for n = 0,1,2,... (with the nnderstanding that b_i{9,u)) = 0). Eqnation fl27|) is a time- 
independent generalized version of the Brinkman’s hierarchy [22l[28l[29] . The hierarchy was 
introdnced to stndy the approach to a stationary state of a system of noninteracting particles 
snbjected to external forces and noise. We remark that in our case we have a system of 
interacting particles, so that the forces are both external (the driving torqnes) and dne to an 
interaction potential. The eqnation for n = 0 recovers the result that bi{9,u) is independent 
of 9. Noting the scaling of the various terms in Eq. fl27l) with m, we expand bn{9,u)) as 

CX) 

bn{9,u) = '^{^/m)^Cn,k{9,UJ), (28) 

k=0 

with bi{9,u) independent of 9 implying that so is ci^k{9,u)) V k. The only constraint on 
bQ{9,u) being d6* bQ{9,u) = 1, we may withont loss of generality choose Co,fc>i(0,a;) = 0. 
This will prove very usefnl for further analysis, as will be shown below. We now use Eq. fl28l) 
in Eq. fl271) and eqnate to zero the coefficient of each power of The term proportional 
to (\/m) gives simply 


nCn,o{9,uj) = 0, (29) 

which implies that Cnfi{9,u) = 0 for n > 0. The coefficient of the term proportional to 
leads to 

^/{n + y/nTa{9, uj)cn-i,k{9, co) + nCn,k+i(9, ca) = 0 

(30) 

8 















t 



\ 


Figure 1. Flow diagram for the evaluation of the expansion coefficients = 

0,1,2,. .., 6 . Starting from the main diagonal, arrows and different colors denote subsequent 
flows (see text). The elements below the main diagonal are all zero. 


for n,k = 0,1,2,... (with C-i^k{0,u) = 0), where a{9,u) = [G( 6 ') — auj]/T. The system of 


equations (l30|) can be solved recursively, as we detail now. 

3.1. Solution of the system of equations 

Let us hrst consider Eq. for n = 0; we immediately obtain that c\^k{9,uj) is independent 


of 6 for each /c, as we had already inferred above. To proceed, we consider Eq. fl30|) for 


/c = 0 and n = 2, 3,...; since CnfliO, cn) = 0 for n > 0, we get that Cn,i{0, cn) = 0 for n > 1. 
Considering next the equation for fc = 1 and n = 3,4,..., we get Cn^2{9,uj) = 0 for n > 2; 
for /c = 2 and n = 4, 5,... we get 3 (^, 0 ;) = 0 for n > 3, and so on. We thus arrive at 
the general result that Cn,k{9,u)) = 0 W k < n. In Fig. [H we display the coefficients Cn,k in 
a matrix. The result just obtained shows that the matrix is upper triangular. We are thus 
left to consider Eq. (IHUil for n = 1, 2,... and fc > n — 1, or, equivalently, for k = 0,1,2,... 
and ? 7 , = l,2,...,fc + l. In what follows, we will obtain the elements of the main diagonal, 
Cn,n{9,u), then the elements of the hrst upper diagonal, Cn,n+i{9, 00 ), the elements of the 
second upper diagonal, Cn^n+2{9,u)), and so on. Thus, let us begin by studying the case 
n = 1 and k = 0. We have 



In this equation, C2fi{9,u) = 0, while is independent of 9. We thus have a hrst-order 

differential equation for co,o(^, <^), with an unknown constant. The requirement of periodicity, 
i.e., CQfi{9,uj) = Co, 0 ( 6 * + 27r, cu), hxes the value of this constant, and we obtain 



(32) 


( 33 ) 
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where g{0,uj) = d0'a{0' ,u}), and co,o(0,a;) is to be fixed at the end by the normalization 
of bo{9,u)). Having determined Cofi{9,u) and Ci^i(a;), we then obtain recursively the main 
diagonal elements, by considering Eq. (j30|) for n = 2, 3,... and k = n — 1; this gives 


/-^dCn-l^n-l{9,Uj) n - -—-^dCn+l,n-l{9, uj) 

- ag -+ V(« + I)r- ^ - 

+\/nTa(0,w)c„_i^„_i(0,w) +nCn^„(0,uj) ^ 0. 


(34) 


Since c„+i^„_i(6', cn) = 0, we have 


C-n,n{9 , Cj) 


l^n—1 {9 ^ Cj) 

F9 


+ a(0, u)Cn-l^n-l{9■, Oj) 


(35) 


for n = 2,3,.... In particular, for n = 2 the first term within the square brackets is absent, 
since is independent of 9. We note that all the functions Cn^n{9,u)) are proportional 

to co,o(0,a;). 

We now determine the elements of the first upper diagonal. We consider Eq. fl30l) for 
n = 1 and k = 1: 

^ +<=1,2(^) = 0. (36) 

This equation has exactly the same structure as Eq. ea, since C2,i{9,u) = 0, and ci^ 2 {^) is a 
constant independent of 9. At this point, we use the fact that Co,A;(0,a;) = 0 for A; > 1. Then, 
the solution of Eq. fl36|) is simply co,i(6*,a;) = ci, 2 (<^) = 0. Next, by considering Eq. fl30|) for 
n = 2,3,... and k = n, and proceeding similarly, we obtain that all the functions Cn,n+i{9, u), 
i.e., the elements of the first upper diagonal of Fig. [H vanish. 

Next, we determine the elements of the second upper diagonal, beginning by considering 
Eq. (I3U]) for n = 1 and k = 2: 


In this equation, C2,2{9, oo) is known from Eq. (jSS]). Then, from the requirement of periodicity 
of co^2{9,u)), and using cq, 2 ( 0 , 0 ;) = 0, we obtain the solutions 

r 27 T ifl/ dc2,2(B',ui) ci{e\ui) fO 

co, 2 ( 9 ,u) = V2^ — I d 9 'e^^^'’^^ 


/o d0'e9(^'--) 


d9' 


€ 1 , 3 ( 0 ;) = - 


r27r jnidc2M9\j^ g{g',uj) 

Jo dd' ^ 


/y dfl'esf"'-") 


(38) 

(39) 


Again, these functions are proportional to co,o(0,o;). Having determined co ,2 and ci, 3 , we 
obtain recursively the elements of the second upper diagonal, i.e., the functions Cn^n+ 2 , from 


10 

















Eq. (I30|) by considering n = 2,3,... and k = n + 1: 

/—p=dCn-l^n+l{(^, n -.^rpdCn+l,n+l{0,Oj) 

+ VvTa{e, u)Cn-l^n+l{0, + nCn,n+ 2 { 0 , Ul) = 0 . 

With the main diagonal elements already determined, this gives 


( 40 ) 


Cn,n+2(^5^) 


dCn.— 


n—1,72+1 


( 0 , 0 ;) 


89 


+ a( 0 , a;)cn-i,n+i( 0 , 


\/{n + 1)T 9cn+i_„+i(6', cu) 


(41) 


n do 

for n = 2,3,.... In particular, for n = 2, the hrst term within the square brackets is absent 
as 01 ^ 3 ( 0 ;) is independent of 9. Also these functions are proportional to Co,o( 0 ,a;). 

Next, we show that the elements of the third upper diagonal vanish. Considering now 
Eq. (ISU]) for n = 1 and k = 3, we have 

+ Vfa(e, ..Me, = 0, ( 42 ) 

In this equation, 02,3 has been previously determined to be vanishing identically, so that the 
solution of the last equation is simply 00 , 3 ( 0 , 0 ;) = 01 , 4 ( 0 ;) = 0. Then, considering Eq. flSU]) 
for n = 2,3,... and k = n + 2, we hnd that all the elements of the third upper diagonal, 
On,n+ 3 , vanish. 

At this point, the procedure of determining the coefficients o„,fc’s should be clear. All 
the elements of the upper diagonals of odd order vanish, this being equivalent to the fact 
that in the portion of each row above the main diagonal, one element every two vanishes, 
i.e., c„,„+i+ 2 fc = 0 for n,k = 0,1,2 ,.... All the nonvanishing elements are proportional 
to co,o(0,a;). The expressions for the main diagonal elements have already been given in 
Eqs. ([32]), fl3^ and fl35|) . On the basis of the analysis above, we can write down the general 
expressions for the nonvanishing non-diagonal elements as 


Co,2k{0,Uj) — V 2 


r27T dc2,2k(9',‘^) 

Jo dd' ^ 




-9{d,ui) / (\Q' 




dO' 


Cl,l+2k{^) — — 


r27r ip/ dc2,2k{9', 1 ^) a(e'.uj) 
Jo 00' ^ 




IQ \ _ IQ \ I \ V^^C3,i+2fc(0, 1^) 

C2,2+2fc(0,<^) — — y —(^{0,U)Ci^i+2k{^) - - --: 


(43) 

(44) 

(45) 


Cn,n+ 2 fc( 0 : 


— 1 ,72— 1 + 2 /c ( 0 ) 


+ a(0, Uj)Cn-l^n-l-\-2k{9Uj) 


+(n-h 1)T dCn+l,n -l+2k{0, Uj) 


n 


do 


n > 3, 


(46) 
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with A; = 1 , 2 ,.... 

We show schematically in Fig. [T]the flow of the solution up to n = fc = 6 , while that 
for higher values proceeds analogously. As shown, the system fl30|) computes progressively 
each element of the main diagonal, and then the elements of the second upper diagonal, each 
one determined by the knowledge of two previously determined elements, and so on. Each 
element of the matrix is proportional to cq o( 0 , cu), which is fixed by the normalization of 

2 fe( 6 *, cu) = 1. The values of Rf have to be determined self-consistently. 

We end the section by pointing out that for a = 0, or equivalently for cu = 0, the 
equilibrium Gibbs-Boltzmann distribution is simply recovered in our procedure. In fact, 
it is immediate to see that for w = 0 the solution of is 6o(^)0) ~ exp[—J d9G{9)] 
and bn{9,0) = 0 for n > 0. Coherently with this, the solution of fl30l) for cu = 0 is 
Co,o(^)0) ~ exp[—J d 6 *G( 6 *)], with all the other Cn^k{9,0) = 0 vanishing. This can be readily 
obtained from Eqs. (j32i) - (l46ll by the fact that g{27r, 0) = 0. 

4. Illustration for a representative example: The Kuramoto interaction 
potential 

In order to illustrate an implementation of our method, we now apply it to the Kuramoto 
potential. In this case, as noted above, only the first Fourier term with s = 1 needs to 
be taken into account in the interaction potential. For illustrative purpose, let us choose 
a representative Q{uj), namely, a Gaussian: Q{uj) = l/(\/^) exp(—a;^/2), and study two 
physically relevant quantities in the synchronized phase. One is the marginal ^-distribution, 
n{9) = fZ^duj ^(cu) dn f^y’^(9,v,u), i.e., the density profile; using the orthonormality 
of the Hermite functions, one gets 

/ CX) 

duj 0(u)bo(9,u). (47) 

■CO 

The other is the quantity p(9) = fZ^doj Q(uj) fZ^dv v^f^y’^(9,v,uj), which is proportional 
to the local pressure [30] • Using again the orthonormality of the Hermite functions, one has 

p{9) = T [ du Q{u) {y2b2{9,u)+ bQ{9,u'^ . (48) 

We thus need the coefficients bo{9,u)) and b2{9,u), whose evaluation requires truncating the 
expansion (128|) at suitable values fctrunc of k. From Fig. [H we see that knowing C 2 , 2 k allows 
to compute Co, 2 fc, so it is natural to choose the same fctmnc for both b^id^uj) and 62 ( 6 *, w). 

In Fig. [2], we demonstrate an excellent agreement between theory and simulations for 
the density n{9), for given values of {m,T,a). The simulations are performed through 
integration of the 2N equations of motion flT^ by using the algorithm of Ref. [13] and 
timestep 6t = 0.01. From the plots of the figure, it is evident that our analytical approach 
works very well for both small and large values of m. The agreement is confirmed also in 
the plot of the quantity p{9), proportional to the local pressure, as shown in the left panel of 
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Figure 2. Density n{9) in the dynamics (IT^ with ui = 1, Us>i = 0, a Gaussian t/(w), for 
m = 0.25, T = 0.25, a = 0.295, /ctrunc = 12 (left panel), and for m = 5.0, T = 0.25, a = 0.2, 
fctrunc = 2 (right panel). Simulations (points) are for N = 10®; the theoretical predictions 
are denoted by lines. 
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Figure 3. The left panel shows the pressure p{9) for the same parameters as in the left panel 
of Fig. [2] Simulations (points) are for N = 10®; the theoretical predictions are denoted by 
lines. The right panel shows the local temperature T{9) = p{9)/n{9) and its anticorrelation 
with the density n{9). 


a 

0.0 

0.05 

0.1 

0.15 

0.2 

0.25 

Rf (Theory) 

0.829 

0.825 

0.813 

0.789 

0.75 

0.686 

Rf (Simulations) 

0.829 

0.825 

0.812 

0.787 

0.747 

0.680 


Table 1. vs. a obtained in theory and simulations in the dynamics (HU, with 

Ml = 1, Ms>i = 0, a Gaussian t/(a;), for several values of cr at m = 0.25, T = 0.25, fctrunc = 12. 

Fig. [3l where the parameters are the same of those in the left panel of Fig. |2l As a further 
demonstration of the validity of our method, we list in Table [T] the value of Rf obtained in 
theory and simulations {N = 10®) for several a’s and m = 0.25, T = 0.25; again, we observe 
a very good agreement, within numerical accuracies. 
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Figure 4. Density n{9) in the dynamics (IT^ with ui = 0.3, U 2 = 0.7, Us >2 = 0, a Gaussian 
C/(w), m = 0.25, T = 0.25, a = 0.295, fctrunc = 4. Simulations (points) are for N = 10®; the 
theoretical predictions are denoted by lines. 


The ratio p{6)/n{6) gives the temperature T{6). In equilibrium, one has a spatially 
uniform temperature prohle, i.e., T{9) equals the temperature T of the heat bath, 
independent of 6. Then, the spatially non-uniform temperature prohle in the right panel 
of Fig. [3] (where we show the theoretical computation), is a further demonstration that the 
synchronized state is a NESS. The panel also depicts a density-temperature anticorrelation, 
i.e., the temperature peaks at a 6^ at which the density is minimum, and vice versa. This 
phenomenon of temperature inversion occurs in inhomogeneous plasmas (e.g., the Solar 
corona EB. interstellar molecular clouds 1321 ), and is argued mainly by simulations to be 
a generic feature of long-range interacting systems in NESSs [33l[3l]; here, we provide an 
analytic demonstration of the phenomenon. 

To illustrate the generality of our method, we show in Fig. |4] the results of adding a 
cos20 interaction to the Kuramoto potential (ui = 0.3, M 2 = 0.7, Ms >2 = 0). Also in this case 
we see a perfect agreement of the theory with simulations. 

We now discuss the behavior of the truncation order /ctmnc in computing the density 
n{9) as a function of m at a given representative (cr, T); in particular, we point out that Eq. 
(|28|) is an asymptotic expansion in the inverse damping coefficient y/m. Let us consider, 
e.g., parameter values used in this work, i.e., a = 0.295, T = 0.25 and m = 0.25. For 
these values /ctrunc = 18 gives a perfect match with simulation results, as in the left panels of 
Fig. |2]and Fig. |3l The match is no more perfect and gets worse and worse on successively 
including more higher order terms in the expansion; in this case, using the Borel method [23] 
of summing a divergent series circumvents the problem, allowing to correctly compute n{9) 
for a truncation order that in principle could be arbitrarily large. This is expected of an 
asymptotic expansion, and makes us conclude that Eq. fl28l) is an asymptotic expansion in 
y/m. In the Appendix, we give some details on the Borel summation method. 
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5. Conclusions 


We have proposed a novel method to compute the inhomogeneous NESS distribution of a 
wide class of mean-held systems of rotators subject to quenched disordered external drive 
and dissipation. We have demonstrated an excellent agreement between simulations and 
theory for the noisy inertial Kuramoto model of spontaneous collective synchronization, and 
for its two mode generalization. 

Our method is based on a series expansion of the stationary distribution function 
n, cu). First, the velocity dependence of the distribution is separated from the {9,u)) 
dependence by an expansion in Hermite functions of the velocity, with coefficients functions 
of 9 and u (Eq. (j23l) ): in turn, the latter functions are expanded in powers of the inverse 
friction constant (Eq. (l28|) h The second expansion is asymptotic, but we have shown 
that, as is generally the case with asymptotic series, we get the “right” sum by truncating 
at an appropriate order. Furthermore, as mentioned above and as detailed in the Appendix, 
one can apply the Borel summation method to sum the expansion, a method that often sums 
correctly an asymptotic series. We stress that the appropriate order of truncation may be 
found even without resorting to a comparison with simulation data, since computation with 
a larger order of truncation leads to numerical instabilities in the form of oscillations in the 
distribution, as shown in the Appendix. 

We note that our method does not determine if the computed inhomogeneous stationary 
solution is stable for given values of the parameters. For this, it would be necessary to 
perform a stability analysis, which for inhomogeneous solutions is much more complicated 
than that for homogeneous solutions. However, by Ending the inhomogeneous solutions, 
one can theoretically determine the hystheresis loops associated with the presence of non- 
equilibrium first-order phase transitions in the class of models we considered. In fact, the 
knowledge of the stability of the incoherent 6'-independent solution (]22|) as a function of the 
parameters [19], and the determination of the synchronized coherent solution, together allow 
to localize the hysteresis loops in the parameter space. 

To sign off, we want to stress that the method can be applied also to classes of models 
that generalize the one considered in this work [351136] . 
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Appendix: Convergence properties of the density expansion 


Consider an asymptotic power series in the real variable x, 

OO 

A{x) = ^akx'^, (A.l) 

k=0 

and define the partial sum 

n 

An{x) = ^akX^. (A. 2 ) 

k=0 

Being asymptotic means that at any given x 7 ^ 0, one has A„(x) —)■ 00 as n —)■ 00 . In this 
case, one might resort to the Borel summation method by defining the Borel transform of 
A{x) as [23] 

CX> 

^ E (A-3) 

k=0 

If BA{t) converges for any positive t, or, if it converges for sufficiently small t to an analytic 
function that can be analytically continued to all t > 0 , and if the integral 

dt exp{—t)BA{tx) (A.4) 

exists and equals Ab{x) (where the subscript B stands for Borel), then we say that the 
Borel sum of the series on the right hand side of Eq. flA.ip is Ab{x). It is not difficult to see 
that if the original series converges, i.e., if hm„^oo A„(a;) = A{x) < 00 , then Ab{x) = A{x). 
Applying the above formalism to Eq. fl28l) . we get 



boB{.0,uj)= f dt exp(-t) 

k=0 

^ ^ i / /—^ '^0,A;(6') k 


(A.5) 


The last integral is to be computed numerically. One has to truncate the series at a certain 
order k = fctrunc, and to extend the integral over ?/ up to a given value yu, which is chosen such 
that the integrand becomes negligible for y > yM- However, contrary to what happens in the 
original series, we found that the sum in the last integral converges, at least for all |/-values 
smaller than yM that are necessary to compute the integral. We do not know the function to 
which our Borel transform converges, and the corresponding radius of convergence, but the 
numerical results show that our series is Borel summable. The left panel of Fig. I All shows 
the result of computing the density 

/ OO 

du Q{u)boB{d,^^) (A. 6 ) 

■OO 
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Figure Al. Density n{9) in the dynamics (I12|) . with ui = 1, Ms>i = 0, a Gaussian 
m = 0.25, T = 0.25, tr = 0.295. The left panel involves theoretical predictions using the 
Borel summation method with fctrunc = 38, while the right panel involves those using direct 
summation with fctrunc = 22. 


m 

0.0625 

0.125 

0.25 

0.5 

1.0 

k 

^max 

60 

32 

18 

10 

6 


Table Al. For the dynamics, Eq. (ED, with ui = 1, Us>i = 0, and a Gaussian G{oj), the 
table shows the maximum truncation order /cmax in the computation of the density n{6) as a 
function of m at a given representative (cr, T) = (0.295,0.25) for which one observes a perfect 
agreement of the density n{6) in theory and simulations, as in Fig. ([ 2 ). The agreement gets 
worse and worse on successively increasing truncation order beyond fcmax- 


for the same conditions as in the left panel of Fig. [2l by truncating the sum in Eq. flA.5ll 
at fctrunc = 38; the plot coincides with the one shown in the left panel of Fig. [2l On the 
other hand, summing the series fl28ll for n = 0 without resorting to the Borel summation 
method, and then computing the density n{6), the result in the right panel of Fig. lAll shows 
that already for truncation order /ctmnc = 22 of the series, one observes instabilities that 
get worse and worse with further increase of the truncation order (see Table lAll listing the 
truncation order fcmax as a function of m, for the same representative {a, T) = (0.295, 0.25) as 
in Fig. ([2|), up to which one observes a perfect agreement of the density n{9) between theory 
and simulations). We conclude from this analysis that the series (|28l) . although asymptotic, 
is effectively summable by the Borel summation method. 
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